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(-H , Abstract. In this work we propose a new numerical approach to distinguish 

between regular and chaotic orbits in Hamiltonian systems, based on the simultaneous 
integration of both the orbit and the deviation vectors using a symplectic scheme, 
hereby called global symplectic integrator. In particular, the proposed method allows 

OA ' us to recover the correct orbits character with very large integration time steps, small 

energy losses and short CPU times. To illustrate the numerical performances of the 
global symplectic integrator we will apply it to two well-known and widely studied 

^O ' problems: the Henon-Heiles model and the restricted three-body problem. 
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1. Introduction 

Hamiltonian systems exhibit phase spaces where regular and chaotic orbits do coexist. 
This fact makes the problem of the numerical characterization of regular and chaotic 
motion an hard task, notably in systems with many degrees of freedom. For this reason 
scientists developed fast and accurate tools to obtain information about the chaotic 
versus regular nature of the orbits of such systems, and efficiently characterize large 
domains in their phase space as ordered or (weakly) chaotic. 

These methods can be roughly divided into two major groups: Lyapunov-like 
methods, i.e. methods based on the study of the evolution of deviation vectors for 
a given orbit; the methods FLI [1], MEGNO [2], SALI [3] and GALI |1] belong to this 
first class. And Fourier-like methods based on the determination of the frequencies of 
the spectrum of some observable related to a given orbit, for instance FMA [5l E] or the 
Spectral method |7]. In the present work we will be interested in the former class, in 
particular to the SALI chaos detection technique. 

To numerically compute such an indicator, one needs to have a good orbit 
determination, but also to integrate the evolution of two deviation vectors. It is well- 
known that symplectic integration schemes outperform the non-symplectic ones, when 
compared using the same order of accuracy and the same integration step size. This 
is mainly due to the very good energy conservation properties, but also to other first 
integrals. As a consequence, larger step sizes are allowed while still keeping a reasonably 
small energy loss. Thus they allow us to enlarge considerably the time span of the 
numerical simulation, without degradating the goodness of the numerical results, that 
is why symplectic integrators turn out to be essential for the long-term evolution of a 
dynamical system. 

In the same way, chaos detection techniques could also benefit from the use of 
symplectic integrators. The aim of this paper is to show that the symplectic integration 
of the deviation vectors can improve the ability of the previous chaos indicators in 
the characterization of regular and chaotic orbits by correctly identifying the orbit 
behavior, using a larger integration step size. Because our method propose to integrate 
simultaneously both the orbit and the deviation vectors, using a symplectic scheme, we 
hereby name this method global symplectic integrator. 

More precisely, we will accurately determine, using a symplectic scheme and a small 
enough integration time step r, the character, i.e. the regular or chaotic behavior, of 
a large number, Nfot, of orbits using one of the aforementioned chaos detectors. Then, 
all the previous orbits will be reanalyzed using both symplectic and non-symplectic 
integrators schemes, but with larger and larger time steps. We will show that the use 
of the global symplectic integrator will allow us to recover a large percentage of orbits' 
characters with very large time steps compared to the non-symplectic one and using 
small amount of CPU time. 

Our symplectic method will be applied to the classical Henon-Heiles model, and 
also to a well-known problem of celestial mechanics: the restricted three-body problem. 
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where the possibihty of using large integration step sizes is essential for the study of 
secular resonances. 

The paper is organized as follows. In Section 121 we describe our method for 
integrating symplectically the deviation vectors used in the chaos detection techniques. 
For the sake of completeness, we present, in Section |3l a brief introduction to the SALI 
chaos indicator. The results of our method applied to the Henon-Heiles system are 
presented in Section 14. ![ while Section 14.21 will be devoted to present an application to 
the restricted three-body problem. Finally in Section [5] we will summarize our findings 
and we will conclude. 

2. The method 

Let us consider a generic Hamiltonian vector field 

f = JV,H{x) , (1) 

where x = {p^q) & B?"' is the momentum-position vector in the phase space, H{x) is 

the Hamilton function describing the system and ^ = " " j is the standard 

constant symplectic matrix. The solution of ([1]) with initial datum xq can be formally 
written as 

y,(t) = e*^«xo, (2) 

where Lh is the Lie operator, i.e. for all smooth function defined in the phase space we 
have Luf = {H, f}, being {■, ■} the Poisson bracket. 

Given a deviation vector v, its time evolution is described by the tangent map: 

f^=M{^{t))v, (3) 

where M is the Jacobian of the Hamiltonian vector field ([T]), namely M = JV^H, 
evaluated on the solution (p{t). 

Because V^if is a symmetric matrix, the matrix M is an Hamiltonian one, i.e. 
M^J + JM = 0. Hence also the vector field ([3]) is Hamiltonian, with Hamiltonian 
function 

K{v) = ^v^Sv, (4) 

where S = WlH. 

Let us now assume that the function H can be decomposed into the sum of two 
parts each one separately integrable, for instance each one depending only on one group 
of variables 

H{x) = A{p) + B{q). (5) 

Let us observe that in the case of the restricted three-body problem (see Section H^ . 
we will prefer to use a different splitting of the Hamiltonian function: actually the 
function A will correspond to two non-interacting two-body problems that in rectangular 
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coordinates will thus result a function of positions and momenta; while the remainder 
will be by definition the B function. 

Under the above assumptions, a family of suitable symmetric symplectic schemes 
has been presented in [8], hereby called SAB An. Such methods develop the solution 
given by ([2]) using the Baker-Campbell-Hausdorff formula [9] , allowing to approximate 
the true solution by a finite number of compositions of integrable symplectic maps, being 
the error of the approximation a function of the integration step size r. 

More precisely one can find coefficients (cj)j=i_..._n+i and ((ii)j=i,...,n such that the 
map 

is symmetric under the transformation t h-?- —t, is symplectic and is of order 0{t'^^), 
namely there exists an Hamiltonian function H whose exact flow is given by ([6]) and 
moreover H = H + C(r^"). 

Rewriting the deviation vector as -y = {vp,Vq), i.e. identifying its components with 
the natural splitting of the phase space, we can explicitely write the Hamiltonian K as 
follows: 

-v^A{p)Vp + -'og 
where 



Kiv) = -^Aim + l^^Mq)v, , (7) 



A{p)=VlA{p) and B{q) = VlB{q). (8) 

Hence the symmetric symplectic methods SAB An can be used also to integrate the 
evolution of the deviation vectors. 

The aim of this paper is to show that the use of a symplectic integrator to solve 
both the orbit evolution and the tangent equation, i.e. to get the time evolution of the 
deviation vector, can improve the capability of some widely used chaos indicators (such 
as SALI) to determine the character of the orbits in an Hamiltonian system. 

Our claim will be numerically proved for the Henon-Heiles system, using SABA^ 
and SABA2 symplectic schemes in comparison with the non-symplectic 4*'' order 
Runge-Kutta method (see Section H?T|) : while SABAiq and Bulirsch-Stoer will be used 
for the study of the restricted three-body problem (see Section 14. 2p . In the rest of 
the paper the easily implementable SALI chaos detection technique will be used to 
numerically determine the orbit character. For the sake of completeness, the next section 
will be devoted to a brief introduction of the above mentioned chaos indicator. 

3. The SALI chaos indicator 

The Smaller ALignment Index, SALI [3], has been proved to be an efficient and simple 
method to determine the regular or chaotic nature of orbits in conservative dynamical 
systems. Thanks to its properties it has been already successfully applied to distinguish 
between regular and chaotic motion both in symplectic maps and Hamiltonian flows 

(e.g. [mini [El [13]). 
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For the sake of completeness let us briefly recall the definition of the SALI and its 
behavior for regular and chaotic orbits, restricting our attention to symplectic flows. The 
interested reader can consult [Sj [Tl] to have a more detailed description of the method. 
To compute the SALI of a given orbit, one has to follow the time evolution of the orbit 
itself and also of two linearly independent unitary deviation vectors '&i(0),'02(0). The 
evolution of an orbit is given by ([2]), while the evolution of each deviation vector is given 
by the tangent map ([3]). 

Then, according to pj the SALI for the given orbit is defined by 

SALl(t) = mm{\\vi{t) + V2{t)\\ , \\vi{t) - V2{t)\\} , (9) 

where || ■ || denotes the usual Euclidean norm and Vi{t) = infr^, i = 1, 2, are normalized 
vectors. 

In the case of chaotic orbits, the deviation vectors -Oi, V2 eventually become aligned 
in the direction defined by the maximal Lyapunov characteristic exponent [15] (LCE), 
and SALl(t) falls exponentially to zero. An analytical study of SALI's behavior for 
chaotic orbits was carried out in [Tl] where it was shown that 

SALl(t) oc e-^^'-^'^^ (10) 

with ai, cr2 being the two largest LCEs. 

In the case of regular motion, on the other hand, the orbit lays on a torus and the 
vectors vi, V2 eventually fall on its tangent space, following a t~^ time evolution, having 
in general different directions. In this case, the SALI oscillates about non zero values 
(for more details see [ID]). This behavior is due to the fact that for regular orbits the 
norm of a deviation vector increases linearly in time. Thus, the normalization procedure 
brings about a decrease of the magnitude of the coordinates perpendicular to the torus 
at a rate proportional to t~^ and so -Oi, V2 eventually fall on the tangent space of the 
torus. 

The simplicity of SALI's definition, its completely different behavior for regular and 
chaotic orbits and its rapid convergence to zero in the case of chaotic motion are the main 
advantages that make SALI an ideal chaos detection tool. Recently a generalization of 
the SALI, the so-called Generalized Alignment Index, GALI, has been introduced [^IT6]. 
which uses information of more than two deviation vectors from the reference orbit. 
Since the advantages of GALI over SALI become relevant in the case of multi-dimensional 
systems, for the aim of the present paper we will restrict our attention to the SALI. 

The scheme we used to numerically compute the SALI consists in integrating both 
the orbit and the deviation vector using the symmetric symplectic SABA method and 
then compute the indicator according to the definition ([9]). 

4. Results 

In this section, we will present the results of the application of the global symplectic 
integrator to study the orbit behavior of two well-known problems: the Henon-Heiles 
system and the restricted three-body problem. The analysis of the former widely studied 
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Figure 1. Hcnon-Heiles phase space (section x = 0). The energy has been fixed to 
E = 1/8 and three characteristic orbits have been computed on this energy level: a 
regular orbit (black) with initial data a;(0) = 0, j/(0) = 0.55, Px ^ 0.2417 and Py = 0, 
a chaotic orbit strongly diffusing (red) with initial data x{0) — 0, 2/(0) — —0.016, 
Px ^ 0.49974 and Py = 0, and a chaotic orbit slowly diffusing (blue) with initial data 
x{0) = 0, 2/(0) = -0.01344, px - 0.49982 and Py = 0. 



dynamical system, will point out that chaos indicator can benefit from our method, 
regarding the characterization of regular and chaotic motion (see Section I4.ip . The 
second application of our method consists in a classical problem of celestial mechanics, 
namely the restricted three-body problem, that will be presented in Section 14.21 

4.1. Henon-Heiles system 

The Henon-Heiles Hamiltonian system [17] is a well-known and studied model described 
by the following Hamilton function 

H{px,Py,x,y) = - {pI + pI + x"^ + y^) + x^y - -y^ . (11) 

In the phase space, regular and chaotic orbits coexist (see Figure [1]), whose character 
can be accurately and rapidly determined using SALI and integrating both the orbit 
and the deviation vectors using the SABA^ scheme with a time step size r = 10~^. 
The choice of a 4*'' order integrator is motivated by a balance between reliability and 
computation time. Indeed, the smaller the used time step or the higher the integrator 
order, the better the energy preservation, but the longer the computation time. In this 
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Figure 2. Characterization of tlie tliree orbits of Figure [T] a regular orbit (black), a 
chaotic orbit strongly diffusing (red) and a chaotic orbit slowly diffusing (blue). The 
integration step size has been fixed to r = 10"'' and both the orbit and the deviation 
vectors have been numerically integrated using the SABA/^ method. 
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Figure 3. Relative energy loss, maxo<f<T/„ \AE{t)/E{t)\, for SABA4, SABA2 
and RKA integrators as a function of the time steps. Both quantities are given in 
logarithmic scale. Several integration times are selected, T/i„ = {10^,10^,10^}, as 
reported in the legend. 



respect, we decided to adopt a SABA4 integrator with a sufficiently small time step of 
r = 10-1 

In Figure Ej we report the results of the numerical computation of the SALI for 
the three generic orbits of Figure [H Since the characteristic period of the orbits on 
this energy level, is of the order of 10 time units, the integration time span has been 
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Figure 4. Comparison of integration schemes on the regular orbit x — Q^ y — 0.55, 
Px ~ 0.2417 (top panels) and Py — and on the chaotic orbit x — 0, y ~ —0.016, 
Px ~ 0.49974 and py = (bottom panels). Left panels, SALI computed using non- 
symplectic RKA for time steps r e {0.001,0.01,0.05}. Right panels, SALI computed 
using SABAi for the integration of both the orbit and the deviation vectors, for 
the same time steps. Both methods determine correctly the character of the orbits 
independently from the small time step used. 



fixed to 10^ time units. We can observe that the three possible dynamical behaviors, 
regular orbits, chaotic orbits strongly diffusing and chaotic orbits slowly diffusing, are 
well identified. Let us in fact observe that the strong chaotic behavior of the red orbit 
is translated to a quick decrease of SALI to zero. On the contrary, for the black regular 
orbit, SALI remains bounded away from zero. The blue orbit has a particular behavior: 
for quite a long time, up to ~ 30 000 time units, this orbit "follows closely" a periodic, 
thus regular, orbit, and SALI remains positive, but eventually the chaotic character of 
the orbit manifests and the indicator correctly goes to zero. Actually the above orbit is 
close to an unstable periodic orbit. 

To check the robustness of our method with respect to a non-symplectic one, we 
firstly reanalyze the above three orbits using larger time steps. More precisely, we 
numerically computed the SALI indicator using SABA2, SABA^ and RKA integrators 
using time steps for which the energy loss is bounded by sufficient small value, 10~^. This 
can be easily computed using the results of Figure [3l where we report the relative energy 
loss, averaged over 50 randomly chosen orbits, for the integrators SABA^., SABA2 and 
RKA, as a function of the time steps. For sake of completeness, we computed the 
relative energy loss using several integration time spans. Let us emphasize that, the 
energy losses for symplectic integrators are almost the same for all the integration time 
spans for a fixed value of r, on the other hand the non-symplectic scheme RK4 exhibits 
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a strong dependence on the integration time span: the larger the time span is, the worst 
is the energy loss for a fixed value of r. As expected, the larger the time step used, the 
larger the loss of energy. In the following, we fix the largest time steps, r, such that the 
relative energy loss is smaller than 10~^, for an integration time span equal to 10^ time 
units. Under these assumptions we get: Tmax ~ 0.3 for SABA4, Tmax ~ 0.16 for SABA2 
and Tmax ~ 0.08 for RKA. As a result for the Henon-Heiles system, our symplectic 
scheme allows time steps four times larger than a non-symplectic one of the same order, 
as will be shown hereafter. 

With these maximal realistic time steps in mind, we reanalyzed the three orbits of 
Figure [2] with time steps r G {10~^, 0.01, 0.05} for which the energy loss is small enough, 
even for the RKA integrator. All the integration schemes behave almost equally for all 
the small used time steps, as it can observed from the results reported in Figure HI 

The next step of our comparison is to consider a larger portion of phase space to 
capture information on the characterization of the global dynamics using symplectic and 
non-symplectic integration schemes. We hence consider N^eg = 100 randomly chosen 
regular orbits and N^ha = 100 randomly chosen chaotic orbits, whose behavior has been 
accurately determined using a sufficiently small step size r, namely r = 10~^, and a 4*^^ 
order symplectic scheme. Then we compute as a function of the used step size, how 
many orbits are correctly characterized by the RKA non-symplectic integration, and by 
the SABA2 and SABA4^ symplectic integration for both the orbit and the deviation 
vectors. 

On the one hand, the results reported in Figure [5] clearly show that, for regular 
orbits, almost all schemes are able to recover nearly 100% of regular orbits. Concerning 
the step sizes. Figure [5] shows that SABA^ characterizes correctly the dynamics of the 
same percentage of orbits with 4 times larger time steps than RK4, while SABA2 
does the same with 2 times larger time steps than RK4, but observe that its order is 
half the one of RK4. On the other hand, the same conclusion holds for the chaotic 
orbits. Moreover, let us note a general trend of RKA to overestimate the number of 
regular orbits when increasing the time step. Thus, our symplectic scheme is largely 
better than RKA to characterize regular orbits, regarding both the used time step and 
the integration order. Indeed, SABA2 symplectic integrations reveal to be nearly as 
reliable as SABA4 ones, especially for small step sizes. 

In Figure [6] we report the CPU times averaged over fifty randomly chosen orbits, 
Tcpu, as a function of the time steps. One can clearly observe that the symplectic 
schemes are faster than RK4 using the same step size, more precisely SABA4^ is 
almost 2.5 times faster than RKA. Once again this fact illustrates the good numerical 
performance of our method which is less time consuming in itself but also allows larger 
time steps, reducing again the computational time. 

To quantify the efficiency of our symplectic scheme, we introduce the following 
efficiency indicator: e = p\logiQ{\AE/E\)\\\ogiQ{Tcpu)\, whose dependence as a 
function of the time step r is represented in Figure [3 Let us observe that the larger is e, 
the better is the integration scheme. For small integration time steps, the non-symplectic 
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Figure 5. Global comparison between RKA non-symplectic scheme and SABA/^^ 
SABA2 symplectic schemes for a large portion of the Henon-Heiles phase space: 
Nreg — 100 regular orbits and Ncha — 100 chaotic orbits are considered. Top panel: 
Percentage of correctly identified regular orbits for increasing time steps {p = correctly 
identified orbits / total number of orbits). Bottom panel: Same as top panel for chaotic 
orbits. 



RKA scheme dominates, mainly because its energy loss is measly (see Figure [3]). On 
the other hand, for time steps r larger than 0.1, the situation is reversed: the global 
symplectic method dominates, as its energy loss and CPU times are both relatively 
small (see Figure [3] and [6]) , that results in a quite large efficiency. 

As a result, it appears that the use of a symplectic integrator to compute the time 
evolution of the deviation vector can improve the capability of a chaos indicator such as 
SALI to determine the character of the orbits in a Hamiltonian system. Furthermore, 
we showed that our global symplectic integrator allows larger time steps without energy 
loss and saves a considerable amount of computation time. This possibility turns out to 
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Figure 7. The efficiency e — p\ \ogiQ{\AE / E\)\ \ \og^f^{Tcpu)\ of the global symplectic 
method {SABA4 and SABA2) and the non-symplectic one (RKA). The scales are 
logarithmic. See text for detail. 



be essential for the study of real problems, where the time scales are fixed by physical 
constraints, as for instance in the case of secular resonances in the N-body celestial 
mechanics problem, as shown in the next section. 

4-2. Restricted three-body problem 

The problem of three celestial bodies interacting with each other gravitationally is a well- 
known problem in celestial mechanics. In the following we will consider the restricted 
problem. More precisely one of the bodies, hereby called the Sun, is the largest mass 
around which the two other bodies are assumed to evolve; the second body, possessing an 
intermediate mass hereby called Jupiter, orbits around the Sun on a circular orbit, while 
the third body, hereby called the asteroid, has a negligeable mass and is moving on an 
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inner orbit "between "the Sun and Jupiter. Assuming this geometry, Kozai [IB] showed 
that a highly inchned asteroid perturbated by Jupiter is characterized by a coupled 
variation in its eccentricity e and inclination i, in such a way that H = y/a{l — e^) cos(i) 
is a constant, a being the asteroid's semi-major axis. This dynamics is often referred to 
as Kozai resonance (e.g. |19j). 

Indeed, the restricted three-body problem can be reduced to two degrees of freedom 
after short-period averaging and node reduction (see for instance [20]). Assuming that 
the outer giant planet is on a circular orbit, this problem is integrable, and its dynamics 
can be represented on the phase space (ecosw, esinw) where u is the argument of the 
pericentre of the massless body (see [H], [22] for more detail). Such a representation 
is given in Figure [8] where we plotted several trajectories of the small body obtained 
by the numerical integration of the Hamiltonian equations associated to the democratic 
heliocentric formulation (see for instance [23]) of the three-body problem: 
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3 

To realize the representation of Figure [HI we used a 10*'^ order SABA integrator 
with a time step sufficiently small of r = 10^^ years. This phase space corresponds 
to H = 0.41833, which means an inclination of 45° for the circular massless body 
represented on the center of the plot. As one can see, for that large inclination 
value, circular orbit corresponds to an unstable equilibrium point. A separatrix divides 
the phase space in three parts: two regions are characterized by the libration of u 
respectively around 90° and 270° and a third one corresponding to a circulation of this 
angle. The two stable equilibria enclosed by the separatrix are referred to as Kozai 
equilibria. Hence, a massless body initially on a circular orbit will suffer from large 
variations in eccentricity, since its real motion (short periods included) will stay close 
to the separatrix of the reduced problem. 

Let us note that these perturbations on a small body at high inclination are 
secular, which means that they operate on extremely long time scales (about 10, 000 
years) compared to the annual orbital periodic variation of the bodies. Symplectic 
integrators are recommended for that kind of problem, as they possess good energy 
conservation properties and are efficient for large time steps. Our method for 
integrating symplectically the deviation vectors is then particularly suitable for the 
chaos determination of such a problem. 

In Figure |9l we report the results of the numerical computation of the SALT with 
SABAio integrator (r = 10~^ years) for a regular orbit (e = 0.4, i.e. close to one of 
the stable Kozai equilibria) and a chaotic orbit (e = 0.001, i.e. close to the unstable 
equihbrium). Our global symplectic integrator scheme reproduces correctly the expected 
character of both orbits in less than 40, 000 years. 

Concerning the global behavior of the phase space, it can be thoroughly described 
by resorting to a chaos indicator technique along a cross section of the phase space. 
We choose to represent the SALI values as a function of the initial values for esino;, as 
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Figure 8. Restricted tliree-body problem pliase space, corresponding to H = 0.41833, 
and reproducing the Kozai resonance (see the text for further details). Initial conditions 
for the calculation are the following: a = 0.35, agiant = 1, egiant = 0, igiant = where 
the subscript giant refers to the outer giant planet. 




Figure 9. SALI characterization of two orbits of Figure |8] with a symplectic SABAio 
scheme: a regular orbit at e = 0.4 (red) and a chaotic orbit at e = 0.001 (blue). 



shown in Figure HOl The chaos along the separatrix is clearly visible for values of e sin u 
close to and to 0.5. 

In order to compare our method with a non-symplectic scheme, we tried to 
reproduce these results with the Bulirsch-Stoer integration method. The first comment is 
that this adaptive stepsize method requires computation times up to three times longer 
than the global symplectic integrator. Secondly, it appears that this non-symplectic 
method seems unable to characterize correctly the orbits, regardless of the time step 
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Figure 10. SALI characterization of the global behavior of Figure [8] as a function of 
esincj, after 100,000 years. In order to avoid useless computation time, SALI values 
have been fixed to 10~^^ when reaching this threshold. The influence of the separatrix 
is obvious and well identified by our symplectic SABAiq scheme. 
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Figure 11. Norm of the deviation vector computed with SABAiq symplectic scheme 
(red curves) and Bulirsch-Stoer method (blue curves) for the two orbits of Figure HI 



used. Indeed, even if the orbit is correctly described by the Buhrsch-Stoer method, 
the deviation vectors do not enable us to distinguish a regular behavior from a chaotic 
one. To show this, we reported in Figure [TT] the norm of the same deviation vector 
integrated with our symplectic scheme (red curves) and with the Bulirsch-Stoer method 
(blue curves). Using the SABAiq method, the norm of the deviation vector grows 
very rapidly for the chaotic orbit identified in Figure M On the contrary, no significant 
deviation is observed between the two orbits of Figure [9] for the non-symplectic method. 
Maybe, the integration on a longer time span with the Bulirsch-Stoer method would 
eventually reveal the chaotic behavior of one of these two orbits. Anyway, given that 
the energy is not fully conserved, one could not reliably trust this result, even more on 
long time interval. 
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5. Conclusions 

In this work, we proposed a new metliod for the detection of regular and chaotic orbits 
in Hamiltonian systems, based on the integration of the deviation vectors used in chaos 
detection techniques, using symplectic algorithms. Our method has been tested on 
two well-known models, and the results clearly demonstrate that it outperforms non- 
symplectic ones. 

Concerning the Henon-Heiles system, it appears that, for large time steps, non- 
symplectic integrators tend to detect an excessive number of chaotic orbits, while the 
global symplectic integrator is able to identify correctly the character of nearly all orbits 
for larger time steps, up to four times larger than non-symplectic ones. Moreover, due 
to his symplectic properties, we showed that our method ensures a very small energy 
loss even on very long time spans. Let us emphasize that the possibility to use larger 
time steps saves a considerable amount of computation time. 

This possibility turns out to be essential for the study of the Kozai resonance in the 
restricted three-body problem, where the secular orbital changes operate on extremely 
long time scales. Once again, the influence of the separatrix of this problem is well 
identified by our global symplectic integrator. On the contrary, the Bulirsch-Stoer non- 
symplectic method, even with smaller time steps, seems unable to distinguish between 
regular and chaotic motion, at least on the same integration time span. A possible 
reason for this behavior could be the accumulation of numerical errors introduced by 
the integrator and a significant energy loss, disadvantages which are avoided using our 
symplectic scheme. 

We are confident that our findings should be generic for a large class of Hamiltonian 
systems. Thus we encourage scientists working on chaos indicators to perform symplectic 
integrations of both the orbit and the deviation vectors using the global symplectic 
integrator, as proposed in the present paper, whenever the Hamiltonian is of the form 
H{x) = A['p) +B{q), or generically it can be divided into two parts, each one separately 
integrable. Computation time and reliability of the results could thus benefit a lot from 
this procedure, as we clearly demonstrated above. 
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